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I. 



INTRODUCTION 



A. PROBLEM STATEMENT 

During a military deployment, troops, equipments and 
supplies must be transported from ports of embarkation to 
ports of disembarkation. Generally, the standard modes of 
transportation used in this operation are trucks, trains, 
airplanes and ships. Because of the limited amount of 
available resources and transport assets, planning becomes 
essential for a successful deployment. During peacetime, 
cargo transportation can be routinely scheduled and the 
normal criterion for a deployment plan is its cost (or 
operating expense). However, during a period of conflict 
(or crisis), expenses become secondary and it is more 
important to transport the troops and cargoes to their 
destinations as fast as possible. 

This study restricts itself to the problem of 
constructing an optimal deployment plan which employs only 
sealift assets. Many of today’s naval deployment plans are 
constructed manually and in an ad hoc manner. This process 
is quite time-consuming and does not guarantee to produce 
even a near optimal plan. 
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B . BACKGROUND 



Research in deployment planning for both industrial and 
military applications has been concentrated on constructing 
deployment plans which minimize operating costs. In a 
survey article, Ronen [Ref. l:pp. 119-126] describes the 
various modes of operation for cargo ships and provides a 
classification scheme for ship routing and scheduling 
models. In a more recent article, Brown et al . [Ref. 2:pp. 
335-346] present and solve the crude oil tanker scheduling 
problem formulated as an elastic set-partitioning model. 

On the military side, Goodman [Ref. 3] formulates the 
problem of scheduling the naval surface combatants of the 
Atlantic Fleet as a generalized set-partitioning problem. 
The resulting constraint matrices in both formulations of 
Brown et al. [Ref. 2] and Goodman [Ref. 3] have a large 
number of columns which must be generated beforehand and 
correspond to all feasible ship schedules. In a Naval 
Postgraduate School master's thesis, Collier [Ref. 4] 
formulates the deployment planning problem employing four 
standard modes of transportation as a linear programming 
problem, and solves it by the MPS III Mathematical 
Programming System developed by Ketron Management Science, 
Inc. [Ref. 5]. Related to Collier's study, Lally [Ref. 6] 
uses the General Algebraic Modelling System, GAMS [Ref. 7], 
to solve the problem of minimizing the number of sealift 
assets required to carry out a given deployment plan. 
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C. OBJECTIVE 



In previous formulations of the deployment or ship 
scheduling problem, the primary objective is to minimize 
cost which is the most appropriate for peace-time military 
operations and for industry. This thesis addresses the 
same problem, but with a different objective: to minimize 
the duration of the deployment. In particular, it 
considers the construction of schedules for sealift assets 
to transport cargoes from their ports of origin to their 
ports of destination in the minimum length of time. 
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II. PROBLEM FORMULATION 



To formulate the crisis deployment problem, the 
following data are assumed to be given: 

1. The ports of embarkation and disembarkation for each 
cargo 

2. The distances between ports 

3. The number of ships with their speed 

4. The compatibility between each ship and each cargo 

When a ship is compatible with a cargo, we mean that 
the ship is compatible with both the cargo and its ports of 
embarkation and disembarkation. Therefore, in constructing 
the compatibility information one has to consider, for 
example, the ship and cargo type as well as the ship draft 
and the channel depth of both ports. 

It is assumed that all cargoes are configured into full 
shiploads. This implies that when a ship picks up a given 
cargo, it must deliver it before any other cargo can be 
picked up. Therefore, the ship must travel to the port of 
disembarkation directly from the port of embarkation. 

A. MATHEMATICAL FORMULATION 

The problem of scheduling sealift assets in a crisis 
situation can be formulated as a variation of the standard 
set-partitioning model as follows. 
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Indices : 



i - indexes shiploads of cargoes, where 
i=l,2,...,m and in is the number of 
shiploads ; 

j - indexes ships, where j=l,2,...,n and n 
is the number of ships available; 

k - indexes a ship schedule. 

Data : 

Sj k - a binary vector representing the kth feasible 

schedule for ship j . The i* h component of Sj k , 
denoted by Si j k , is; 

1, if the k t h feasible schedule 

for ship j includes shipload i; 

Si J k = < 

0, otherwise; 

tj k - the completion time of schedule Sj k . 

(The calculation of tj k is described in a 
numerical example below. ) 

Decision variables: 

1, if the k* h feasible schedule 
for ship j is selected for the 
xj k = < deployment; 

0, otherwise. 
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Problem PI; 



min [ max ( Z tik *xu Z tnk*x D k ) ] 

k k 

subject to: 
n 

Z Z Sijk * xj k £ 1 for i=l,...,m 

j=l k 

Z xj k Si for j=l , . . . , n 

k 

Xj k £ (0,11. 

The term Zk tj k *xj k in the objective function 
represents the time for ship j to complete its assigned 
schedule. Therefore, max { Zk ti k *xi k , . . . , Zk tnk*Xok 1 
represents the completion time of the longest schedule in 
the deployment plan. Since the deployment is considered 
completed only when all cargoes are delivered to their 
destinations, the completion time of the longest schedule 
becomes the length of the deployment plan, which is to be 
minimized. The first set of constraints in Problem PI 
ensures that all cargoes are picked up by at least one ship 
and the second guarantees that at most one schedule is 
assigned to each ship. 

In addition, the objective function of Problem PI is a 
nonlinear convex function since it is a point-wise maximum 



( 1 ) 

( 2 ) 
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of a set of linear functions. However, Problem PI can be 
transformed into a linear problem as follows. 

Problem P2 : 

minimize z 

subject to: 
n 

I I Sijk ‘ xjk £1 for i=l,...,m; 

j=l k 

I xjk £ 1 for j=l,...,n; 

k 

z - I tj k * xjk £ 0 for j=l,...,n; 

k 

Xj k E (0,1). 

Note that the last set of constraints defines the objective 
function of Problem PI. 

B . AN EXAMPLE 

To illustrate the above integer programming 
formulation, consider the deployment problem depicted in 
Figure 2.1. There are two ships, Ship 1 and Ship 2, 
available for the deployment. Initially, Ship 1 is docked 
at Port PI and Ship 2 at Port P2. There are 3 shiploads 
whose ports of embarkation and disembarkation are given in 
Table 2.1. Assume that both ships are compatible with all 
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ports and cargoes. The lines connecting various ports in 
Figure 2.1 represent possible movements between ports. The 

TABLE 2.1 

A LIST OF POE AND POD FOR THE SHIPLOADS 



SHIPLOAD 


POE 


POD 


1 


1 


1 


2 


2 


1 


3 


3 


2 



POE ' POD 




Figure 2.1 Data for the Deployment Problem 
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numbers adjacent to the lines represent the travel times 
for both ships, i.e., they have the same speed (assumed 
constant regardless of cargo loading) . 

Consider a schedule for Ship 1 which includes picking 
up cargoes 1 and 3. The binary vector representing this 
schedule has 3 components (since there are 3 shiploads) 
with the first and the third components having the value 
one and the second component the value zero. To carry out 
this schedule, Ship 1 can use one of the two possible 
routes: one which picks up cargo 1 first and then cargo 3 
and the other which picks up the cargoes in the reverse 
order. Using the time given in Figure 2.1, the first route 
requires 8 days to complete and the second requires only 7 
days. Since the objective is to minimize the completion 
time of the longest schedule, the completion time of 7 days 
is assigned to this schedule, i.e., ti 4 =7. In general, 
the completion time tj k is the time for ship j to carry out 
schedule k using the shortest route. Tables 2.2 and 2.3 
display all possible schedules along with their completion 
times for Ships 1 and 2, respectively. Note that the 
schedule discussed above is the schedule Si 4 in Table 2.2. 

The optimal deployment plan for this example consists 
of two schedules: Si 2 for Ship 1 and S 2 4 for Ship 2, and 
requires 7 days to complete. In terms of decision 
variables, xi 2 and X 2 4 equal one and all other variables 
equal zero. 
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The explicit formulation of the above example is given 



in Figure 2.2. 



TABLE 2.2 

POSSIBLE SCHEDULES FOR SHIP 1 





SCHEDULES 


Si i 


Si 2 


Si 3 


Si A 


Si 3 


Si 6 


Si 7 


COMPLETION 

TIME 

(IN DAYS) 


3 


6 


2 


7 


10 


8 


12 


SHIPLOAD 1 


1 


0 


0 


1 


1 


0 


1 


SHIPLOAD 2 


0 


1 


0 


0 


1 


1 


1 


SHIPLOAD 3 


0 


0 


1 


1 


0 


1 


1 



TABLE 2.3 

POSSIBLE SCHEDULES FOR SHIP 2 





SCHEDULES 


Sz i 


S2 2 


S 2 3 


S2 A 


S 2 3 


S 2 6 


S2 7 


COMPLETION 

TIME 

(IN DAYS) 


4 


4 


2 


1 


8 


8 


12 


SHIPLOAD 1 


1 


0 


o- 


1 


1 


0 


1 


SHIPLOAD 2 


0 


1 


0 


0 


1 


1 


1 


SHIPLOAD 3 


0 


0 


1 


1 


0 


1 


1 
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minimize z 
subject to: 

Xll+ Xl4+Xl3 + Xl 7 + X2 1+ X2 4+X2 3 + X2 7 11 

X 12 + X 13 +X 16 +X 17 + X 22 + X 23 +X 26 +X 27 1 1 

X1 3 +X1 4 + X1 6+X1 7 + X2 3 +X2 4 + X2 6+X2 7 1 1 

Xl 1 +Xl 2 +Xl 3 +Xl 4 +Xl 3 +Xl 6 +Xl 7 ^ 1 

X2 1 +X2 2 +X2 3 +X2 4 +X2 3 +X2 6 +X2 7 £ 1 

Z “ [ tl 1 *Xl 1 +tl 2 *Xl 2 +t 1 3 *Xl 3 + tl 4 *Xl 4 +tl 3 *Xl 3+tl 6 *Xl 6 +tl 7 *Xl 7 ] 1 0 
Z ” [ t 2 1 *X2 1 +t2 2 *X2 2 +t2 3 *X2 3 +t 2 4 *X2 4 +t2 3 *X2 3 +t 2 6 *X2 6 +t2 7 *X2 7 ] 1 0 



Xl l , 


Xl 2 , 


Xl 3 , 


Xl 4 , 


Xl 3 , 


Xl 6 , 


Xl 7 


£ 10,1) 


X2 1 , 


X2 2 , 


X2 3 , 


X2 4 , 


X2 3 , 


X2 6 , 


X2 7 


£ 10,1) 



Figure 2.2 A Formulation Of The Example Problem 
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III. A SEQUENTIAL SOLUTION TECHNIQUE 



The solution procedure presented below addresses 
Problem PI (and P2) indirectly. This procedure takes 
advantage of the fact that there exists a simpler problem 
which is related to Problem PI (and P2) . By sequentially 
solving a number of these simpler problems, one can arrive 
at a solution to Problem PI (and P2). 

A. A RELATED PROBLEM 

In certain situations, it is not so critical that the 
planner obtains a deployment with the minimum duration. 
Instead, the duration of the deployment, say r days, has 
been set by the top command and the planner only has to 
find a feasible plan which can be completed within this 
given length of time. To formulate this problem, define: 

Kj (t) = ( k : Sjk is a feasible schedule for ship j 
and t i k < t t , 

That is, Kj (r) is the set of schedules for ship j which 
can be completed within r days. Then, we have the 
following problem which we refer to as the feasibility- 
seeking problem. 
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Problem P3 ( x ) : 



m 

min Z wi 
i=l 

subject to: 
n 

wi + Z ( Z Si j k * xj k ) £ 1 for all i=l,...,m; 

j = l keKj (x) 

Z xj k Si for all j=l,...,n; 
k c Kj ( x ) 

wi , Xj k £ {0,11. 

where wt is an auxiliary variable to indicate whether or 
not shipload i will be left undelivered by the deployment 
plan. If the optimal solution to Problem P 3 ( x ) is greater 
than zero, it means that x is infeasible. In this case, 
one or more shiploads must be left undelivered or 
additional assets are required to obtain a plan which can 
be completed in x days or less. Thus, if x is a feasible 
duration, Problem P3(x) will produce a feasible plan. 

Note that Problem P3(x) is parameterized by x . If the 
minimum duration for a deployment plan, x* , is known, then 
the solution to Problem P3(x*) solves Problem PI (and P2) 
as well. Otherwise, by varying x and resolving Problem 
P3(x) in a systematic manner, one can obtain a solution to 
Problem PI (and P2). A strategy for searching for the 
minimum feasible duration x* is discussed in Chapter IV. 
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To illustrate the feasibility-seeking problem, consider 
the deployment problem presented in Chapter II. Assume 
that the planner is told to construct a plan with a 
duration of 8 days. 

Then , 

Ki ( 8 ) = { 1, 2, 3, 4, 6 ) and 
K 2 ( 8 ) = I 1, 2, 3, 4, 5, 6 I, 

that is, the eligible schedules for this plan with a 
completion time of 8 days or less are those listed Tables 



3 . 1 and 3 . 


2. In this case, 


the 


optimal objective function 


value for 


Problem P3(8) 


is 


zero, because 8 


days is a 


feasible 


duration. Each 


of 


the following 


pairs of 


schedules 


for Ships 1 


and 


2 : ( Si 1 , S 2 6 ) , 


( Si 2 , S 2 4 ) , 



(S13 , S 2 3 ) , ( Si 4 , S 2 2 ) , and (Si&,S2i), constitutes a 

deployment plan that can be completed within 8 days. 

Similarly, if one solves Problem P3(x) with x equal to 
7 days (the optimal duration) , the optimal objective 
function value is still zero, and the pairs (Si2,Sz4) and 
( Si 4 , S 2 2 ) are the only feasible deployment plans. 

B. A COLUMN GENERATION APPROACH TO THE FEASIBILITY-SEEKING 

PROBLEM 

Since the feasibility-seeking problem searches for a 
feasible deployment plan and does not have a real objective 
function to optimize, one expects that the relaxation of 
the integrality restriction would not produce fractional 
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solutions too often. This observation is corroborated by 
the computational result presented in Chapter IV in which 
integer solutions are obtained for over 90 per cent of the 
problems. Henceforth, we treat Problem P3{t) as a linear 
programming problem. 

As a linear program. Problem P3(x) has many columns. 
To avoid generating these columns a priori, we apply the 
column generation technique, i.e., the Dantzig-Wolf e 
decomposition method [see, e.g.. Ref. 8], to Problem P3(x), 
and the following decomposed system is obtained. 



TABLE 3.1 



ELIGIBLE SCHEDULES FOR SHIP 1 WHEN x=8 





SCHEDULES 


Si i 


Si 2 


Si 3 


Si 4 


Sl 6 


COMPLETION 

TIME 

(IN DAYS) 


3 


6 


2 


7 


8 


SHIPLOAD 1 


1 


0 


0 


1 


0 


SHIPLOAD 2 


0 


1 


0 


0 


1 


SHIPLOAD 3 


0 


0 


1 


1 


1 
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TABLE 3.2 



ELIGIBLE SCHEDULES FOR SHIP 2 WHEN x=8 





SCHEDULES 


S2 1 


S 2 2 


S 2 3 


S2 4 


S 2 3 


S2 6 


COMPLETION 

TIME 

(IN DAYS) 


4 


4 


2 


7 


8 


8 


SHIPLOAD 1 


1 


0 


0 


1 


1 


0 


SHIPLOAD 2 


0 


1 


0 


0 


1 


1 


SHIPLOAD 3 


0 


0 


1 


1 


0 


1 



Master Problem (MP) : 



m 

min I wi 
i = l 



subject to: 



n 

2 2 S ) J k * Xj k + Wi £ 1 

j=l keLj (x) 



2 xj k <. I 

keLj (x) 



0 5 xj k £ 1 for 



Subproblem for ship 1 ( SPl(i) ) : 

m 

k' = arg min ( vj + 2 Sijk * 

keKj (x) i=l 
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for 



for 



all 



i = l , 
3=1, 
3, k 



,m (1 



, n (2 



where ui is the dual variable corresponding to constraint 
set (1), i.e., the cargo (shipload) constraints, and vj is 
the dual variable corresponding to the constraint set (2), 
i.e., the ship constraints. We refer to ui as the i* h 
cargo dual and vj as the j th ship dual. 

The column generation technique starts with an initial 
set of feasible schedules, Lj (x), for each ship j. This 
initial set Lj (x) may be an empty set. The master problem 
is solved and the dual variables ut and vj are obtained. 
From this set of cargo and ship duals, one or more 
subproblems are solved thereby generating additional 
schedules (columns), Sj k - , which are subsequently added to 
the set Lj (x) . The master problem is then resolved with 
the additional schedules (columns) and the new cargo and 
ship duals are ained. The cycle then continues until 
the objective junction value of Problem SP(j) is 
nonnegative for all j, i.e., all schedules have nonnegative 
reduced cost. This signifies that optimality is achieved. 
Figure 3.1 illustrates the cycling between the master and 
subproblem in the column generation technique. 

As stated above, the subproblem is unnecessarily hard. 
In theory, it is not necessary to add schedules (columns) 
with the most negative reduced cost to the master problem. 
Any schedules (columns) with negative reduced cost would 
suffice. The following subproblem produces negative 
reduced cost schedules for the master problem. 
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Master Problem 
m 

minimize Z wi 

i=l 

subject to: 
n 

Z Z Sijk * xj k + wi £ 1 , for i=l,...,m 

j=l ktLj { x ) 

Z xj k £ 1, for j=l,...,n 

keLj (t) 

0 £ xj k <. 1 


/ 


\ 




The subproblem 




The master 


produces a new 




problem produces 


column , Sj k > , 




cargo and ship 


for the master 




duals for the 


problem, i . e . , 




subproblem. 


Sj k * is added 


Subproblem 




to Lj (x) . 








\ 


/ 




m 




k ' = arg 


min { vj + Z Sijk * ui ) 


keKj (x) i = l 





Figure 3.1 - The Column Generation Technique. 
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Subproblem SP2(i) : 



For ship j find an index k' such that k' t Kj (t) and 
m 

vj + Z Si J k * Ui <0. 
i=l 

If k' solves Problem SP2(j), k' is an acceptable 
schedule. For details concerning the generation of 

acceptable schedules, the reader is referred to a related 
Master's Thesis by LCdr Svein Buvik [Ref. 9]. 
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IV. IMPLEMENTATION AND COMPUTATIONAL RESULTS 



To implement the column generation procedure we 
modified the revised simplex code described in Ref. 10, to 
solve the master problem. In this modification, we allow 
the algorithm to restart from the last optimal solution 
after one or more new schedules (columns) have been added 
to the master problem. Since the set partitioning problem 
is usually degenerate, we also reinvert the basis at every 
ten iterations. As for the subproblem, we employ the 
algorithm developed by Buvik [Ref. 9]. Both the master and 
subproblem algorithms are written in FORTRAN 77 and 
compiled by the IBM VS FORTRAN compiler. All runs were 
performed on an IBM 3033 AP computer at the W.R. Church 
Computer Center of the Naval Postgraduate School. 

A. PROBLEM DATA 

For our experimentation below, we consider the 
deployment scenario in which cargoes must be moved from the 
ports along the east coast of United States to ports in 
Europe. Table 4.1 lists approximate distances between 
various ports. The number of shiploads for our deployment 
problems are varied between 5 and 50 and the list of all 50 
shiploads along with their POE's and POD's are given in 
Table 4.2. The number of ships assigned to the deployment 
are assumed to be between 2 and 30 ships and the initial 
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location of all 30 ships are given in Table 4.3. The speed 
of all 30 ships are between 18 and 25 knots, and on the 
average a ship is compatible with 75 % of the shiploads. 



TABLE 4.1 

DISTANCES BETWEEN PORTS 
(IN NAUTICAL MILES) 



PORTS 


Ham. 


Wilh. 


Rot . 


Antw . 


Chb . 


n. y. 


4030 


3950 


3790 


3775 


3520 


Norf . 


4340 


4260 


3490 


4075 


3800 


Chari . 


3650 


4560 


5390 


4370 


4090 


Jax . 


4850 


4770 


3470 


4570 


4280 


Pens . 


5390 


3460 


5125 


5110 


4820 



where 



N. Y. 


= 


New York 


Norf . 


= 


Norfolk 


Chari . 


= 


Charleston 


Jax. 




Jacksonville 


Ham. 


= 


Hamburg 


Wilh. 


= 


Wilhemshaven 


Rot . 


= 


Rotterdam 


Antw. 


S 


Antwerpen 


Chb. 


= 


Cherbourg 


Pens . 


= 


Pensacola 
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TABLE 4.2 
LIST OF SHIPLOADS 



SHIPLOAD 


POE 


POD 


1 


1 


1 


2 


2 


1 


3 


1 


3 


4 


2 


2 


5 


3 


2 


6 


3 


1 


7 


1 


4 


8 


3 


4 


9 


1 


5 


10 


2 


5 


11 


2 


1 


12 


4 


1 


13 


2 


1 


14 


2 


2 


15 


3 


2 


16 


2 


3 


17 


2 


2 


18 


2 


4 


19 


2 


3 


20 


2 


5 


21 


2 


3 


22 


3 


2 


23 


3 


1 


24 


3 


2 


25 


3 


4 



POE 

1. New York 

2. Norfolk 

3. Charleston 

4. Jacksonville 

5. Pensacola 



SHIPLOAD 


POE 


POD 


26 


3 


1 


27 


3 


3 


28 


3 


5 


29 


3 


4 


30 


3 


1 


31 


3 


5 


32 


4 


2 


34 


4 


1 


34 


4 


2 


35 


4 


4 


36 


4 


2 


37 


4 


5 


34 


4 


3 


34 


4 


3 


40 


4 


4 


41 


4 


4 


42 


4 


5 


43 


4 


5 


44 


4 


5 


45 


5 


3 


46 


5 


3 


47 


5 


4 


48 


5 


4 


49 


5 


1 


50 


5 


5 



POD 

1 . Hamburg 

2. Wilhelmshaven 

3. Rotterdam 

4. Antwerpen 

5. Cherbourg 
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TABLE 4.3 

INITIAL DISTANCES BETWEEN SHIPS AND PORTS 
(IN NAUTICAL MILES) 



SHIP # 


N.Y. 


Norf . 


Chb . 


Jax. 


Pens . 


1 


154 


245 


550 


720 


1190 


2 


100 


255 


450 


620 


290 


3 


250 


945 


650 


1120 


890 


4 


300 


340 


560 


740 


1890 


5 


250 


320 


990 


2900 


1440 


6 


100 


390 


650 


720 


3290 


7 


245 


120 


300 


475 


975 


8 


245 


230 


300 


475 


975 


9 


200 


•200 


400 


600 


1100 


10 


600 


600 


700 


900 


1400 


11 


150 


100 


400 


575 


1075 


12 


350 


95 


200 


375 


875 


13 


550 


110 


540 


165 


700 


14 


550 


300 


120 


165 


700 


15 


800 


758 


700 


750 


1100 


16 


450 


200 


100 


265 


900 


17 


350 


270 


300 


350 


750 


18 


1240 


1100 


1000 


1100 


1500 


19 


720 


475 


165 


90 


640 


20 


720 


475 


165 


280 


600 


21 


1107 


1500 


900 


905 


1200 


22 


920 


675 


365 


200 


400 


23 


450 


350 


300 


400 


800 


24 


1350 


1250 


1200 


1200 


1600 


25 


1190 


975 


700 


600 


2230 


26 


1190 


975 


700 


600 


2140 


27 


890 


675 


400 


300 


300 


28 


1290 


1075 


800 


700 


300 


29 


900 


700 


500 


400 


600 


30 


1090 


875 


600 


500 


100 
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B. STRATEGIES FOR GENERATING SCHEDULES 



As described in Chapter III, the decomposition process 
iteratively solves the master and subproblem in sequence. 
After having just solved the master problem, the subproblem 
obtains the cargo and ship duals from which it generates 
one or more negative reduced cost columns. At this point, 
there are several possibilities regarding the ship(s) for 
which the subproblem should generate schedules (or 
columns). The first obvious strategy is to generate a 
schedule for Ship 1 in the first iteration, a schedule of 
Ship 2 in the second iteration and so on until a schedule 
for each ship has been generated. At which point, the 
cycle of generating schedules (columns) begins again with 
Ship 1. The second strategy is to generate schedules for 
ship in the descending order of the ship duals, and the 
third strategy is just the reverse of the second strategy, 
i.e., generates schedules in the ascending order of the 
ship duals. The other strategy, which has been considered 
and soon after discarded, generates schedules for all ships 
during each iteration. This strategy tends to generate the 
same schedule for all ships, which seems redundant since no 
two ships can have the same schedule at optimality. In 
fact, there must exists a solution in which no two ships 
are assigned to pick up the same shipload. Based on this 
observation and preliminary experiments, the last strategy 
is discarded. 
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To compare the three strategies discussed above, we 
solved three feasibility-seeking problems at various 
lengths of deployment, x . The first problem has 30 
shiploads and 20 ships, the second problem has 35 shiploads 
and 20 ships, and the last has 40 shiploads and 20 ships. 
In Figure 4.1, we plotted the average cpu times on these 
three problems against the length of deployment. The first 
and third strategies clearly dominate the second. 

C. SOLVING THE MINIMAX PROBLEM 

As mentioned in Chapter III, one can solve the minimax 
problem by sequentially solving the feasibility-seeking 
problem in the following manner. First, pick an initial 
value for x and then solve the feasibility-seeking problem 
at this value x. If the optimal objective function is 
zero, then the value of x is decreased and the feasibility- 
seeking problem is resolved at the new value. Otherwise, 
the optimal objective function is positive, the value for x 
is increased and the feasibility-seeking problem is 
resolved. 

The efficiency of the above algorithm is clearly a 
function of the initial value for x. If the initial value 
for x is close to the optimal, the feasibility-seeking 
problem has to be solved less often. Thus, it is important 
that a good initial value for x is used to start the 
process of increasing and decreasing the value x. 
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CPU TIME 




LENGTH OF DEPLOYMENT 
(IN PERCENT OF THE OPTIMAL LENGTH) 



Figure 4.1 A Comparison of Three Strategies for 

Generating Schedules 
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One lower bound estimate is given by the following 
equation : 

tl = integer part of [(2*ntr - 1) * trmin + itmin] /spmax] 

where 

ntr = average number of cargoes (shiploads) per ship, 
i . e . , 

number of cargoes 

ntr , 

number of ships 

trmin = the minimum travel distance between POD's, 

itmin = the minimum distance between ships’ initial 
positions and POE's, and 

spmax = maximum speed among all ships. 

To understand this bound, note that for each shipload 
assigned to a ship, the ship has to first deliver the 
shipload to its destination and return to pick up the next 
shipload on its schedule. Therefore, the ship has to make 
two trips (or ocean crossings) back and forth between POE's 
and POD's for each shipload, except for the last one for 
which the ship only has to make one trip from a POE to a 
POD. Thus, if ntr shiploads are assigned to one ship, it 
has to make (2*ntr - 1) trips. Since the minimum distance 
between a POE and a POD is trmin, the minimum total 
distance traveled by each ship is (2*ntr 1) * trmin + 
itmin. The first term represents the distance for trips 
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between POE's and POD's and the second term represents the 
distance from ship's initial position to the first POE. 
Then, dividing the total by the maximum speed among the 
ships gives a lower bound for the optimal x. Table 4.4 
displays the value of the lower bound estimate and the 
correspond values of x* , the optimal duration, for 35 
problems. On the average, xl underestimates x* by 40 %. 

If historical data, e.g., data from previous deployment 
exercises, are available, the lower bound estimate xl can 
be improved by using linear regression. For example, using 
the data from Table 4.4, we obtain the following linear 
equation 

test = 15.57 + 0.8 * xl 

where test represent the linear estimate of x* based on xl . 
Figure 4.2 displays the linear estimate of x* graphically. 
Since linear regression minimizes the squared error, some 
lest naturally overestimates x* . Based on Xl and test , we 
implemented the following search algorithm for x* . 

In the algorithm below, the initial estimate, Xi , of 
the optimal duration, x* , is obtained by taking a convex 
combination of the lower bound and the linear regression 
estimates. It is assumed that the convex weight a, is 
chosen so that xi underestimates x* . (Note that this is 
always possible by letting a equals one.) The parameter 5 
equals one time unit which is one day in all our examples. 
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4 

5 

6 

7 

8 

9 

10 

11 

12 

13 

14 

15 

16 

17 

18 

19 

20 

21 

22 

23 

24 

25 

26 

27 

28 

29 

30 

31 

32 

33 

34 

35 



TABLE 4.4 



LOWER BOUND ESTIMATES AND ACTUAL VALUES 
OF OPTIMAL T 



1BER OF 
iPLOADS 


NUMBER 
OF SHIPS 


Tl 


X * 


5 


2 


17 


33 


8 


3 


17 


33 


8 


3 


21 


35 


8 


4 


18 


23 


8 


5 


12 


23 


9 


3 


29 


33 


12 


3 


40 


47 


12 


4 


29 


34 


12 


5 


21 


34 


15 


4 


29 


47 


15 


5 


29 


35 


15 


7 


17 


30 


16 


5 


29 


42 


17 


6 


17 


34 


17 


8 


17 


30 


18 


6 


29 


35 


18 


8 


17 


31 


18 


10 


5 


21 


19 


5 


29 


46 


19 


8 


17 


31 


20 


6 


28 


43 


20 


9 


17 


21 


20 


10 


17 


20 


23 


10 


17 


31 


25 


6 


40 


46 


25 


10 


26 


36 


25 


12 


17 


30 


25 


15 


6 


21 


25 


15 


6 


20 


25 


20 


5 


19 


30 


15 


17 


22 


30 


20 


7 


22 


30 


22 


5 


19 


40 


30 


6 


20 


45 


30 


6 


21 



29 




t l 

Figure 4.2 Lower Pound Estimate versus Actual Values of 



Optimal x. 



Algorithm 



Step 0 : Set Ti = a * tl + (1-a) * Test and 

set k = 1. 

Step 1 : Solve the feasibility-seeking problem, Problem 
P3 ( Tk ) , by the column generation technique. 

Step 2 : If the optimal objective function value equals 

0, stop; tk is optimal. Otherwise, set 
Tk + i = Tk +5 and k = k + 1. Go to step 1. 

In Step 2, the current estimate, Tk , of the optimal 
duration is increased by amount 6. In this manner, the 
current estimate Tk approaches the optimal duration t* from 
below and all of the previously generated schedules remain 
feasible to the feasibility-seeking problems in the 
succeeding iterations. One topic for future research is to 
relax the assumption that Ti must underestimate t* and 
allow Tk to be adjusted in either upward or downward 
direction in Step 2. Table 4.5 summarizes the 

computational results for the above algorithms. In all 
cases, the value for a is 0.7. 
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TABLE 4.5 



COMPUTATIONAL RESULTS FOR THE COLUMN GENERATION TECHNIQUE 



NUMBER OF 
SHIPLOADS 


NUMBER 
OF SHIPS 


RATIO 


Tl 


OPTIMAL 

X 


CPU 

TIME 


8 


5 


1.5 


10 


21 


2.1 


15 


10 


1.5 


10 


21 


9.6 


22 


10 


1.5 


10 


20 


45.5 


30 


20 


1.5 


10 


21 


156.7 


10 


5 


2.0 


19 


21 


2.3 


20 


10 


2.0 


20 


22 


25.0 


30 


15 


2 . 0 


20 


22 


221.3 


40 


20 


2.0 


20 


23 


1216.2 


12 


5 


2.5 


20 


33 


7.9 


25 


10 


2.5 


20 


33 


86.1 
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D. PERCENTAGE OF INTEGER SOLUTIONS 



Since the approach taken in solving the crisis 
deployment problem is the linear relaxation of the minimax 
set-partitioning problem, it is of interest to investigate 
the question concerning the integrality of the obtained 
results. In theory, the linear relaxation of the problem 
does not always produce an integer solution, in which case 
an integer programming algorithm such as the branch and 
bound method must be employed. However, Table 4.6 
demonstrates that an integer programming algorithm is 
necessary for less than seven percent of the problems. 



TABLE 4.6 

NUMBER OF PROBLEMS WHICH TERMINATE WITH OPTIMAL 
INTEGER SOLUTIONS 



T 


# OF PROBLEMS 
SOLVED 


# OF PROBLEMS WITH 
INTEGER SOLUTIONS 


£ 1 . 40*t* 


39 


32 


1 . 30*t* - 1 . 39 * T* 


30 


26 


1 . 20**c* - 1 . 29*t* 


35 


33 


1 . 10 * x* - 1.19*T* 


36 


35 


1.00*x* - 1.09 *t* 


74 


74 


Total 


214 


200 (93.45%) 
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V. CONCLUSIONS AND FUTURE STUDIES 



This study formulates a crisis deployment problem as a 
set-partitioning problem with a minimax objective. An 
algorithm is developed for solving this problem. The idea 
underlying this algorithm is to solve the minimax set- 
partitioning by solving a sequence of simpler, but related, 
feasibility-seeking problems. Each time the feasibility- 
seeking problem produces a better solution for the minimax 
problem. The feasibility-seeking problem is similar in 
form to the minimax problem and both have a large number of 
columns. So to solve the feasibility-seeking problem, the 
column generation technique (as in the Dantzig-Wolf e 
decomposition method) is employed. The computational 
results in Chapter IV verify that this method is effective. 

An important by-product of the above development is 
that the feasibility-seeking problem can also answer the 
question: Can all cargoes be deployed to their final 
destinations in t days? A negative answer to this question 
leads to two natural follow-up questions which provide 
interesting areas for future studies: 

1. How many additional ships are required to deploy all 
cargoes in x days? 

2. If no additional ship is available, which cargoes 
must be left behind? 
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Besides the above areas and the one mentioned in 



Chapter IV, the following areas are also worth studying. 

1. The scenario considered in this study assume that the 
deployment is completed in one phase. In an extended 
period of conflicts, one may want deployment plans in 
several phases (waves) . 

2. Several embellishments to the current model are also 
possible . 

a. Allow the cargoes to arrive at the ports within 
time windows. The current model assumes that all 
cargoes are always available for transport. 

b. Allow cargoes in partial shiploads and in 

compatibility among cargoes, e.g., ammunition 
should not be loaded on same ship with fuel. 

c. Allow for nondeterministic delays in the 

completion times. These delays are due to 

unfavorable weather and/or enemy blockade. 
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APPENDIX 



FORTRAN PROGRAM 



************************************************************ 
* * 



* * 

* * 

* = PROGRAM DEPLOY = * 

* ==================== * 



* 



* 



* * 

* * 

* Date : 23 / 08 / 1988 * 

* * 

* * 



* Key variables 



M 

N 

A 

B 



- number of constraints; * 

- number of variables; * 

- real matrix of dimension M by N containing * 

the coefficients of the M constraints; * 

- real vector of length M containing the right * 

hand sides of the constraints; * 

- real vector of length N containing the * 
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* XB 

* BINV 

* 

* IB 

* 

* U 

* XTIME 

* SEQ 

* ELL 

* K 

* SB 

* SIGB 

* SHIP 

* OBJ 

* MR 

* IT 

* 

* TR 

* 

* SPD 

* TAU 

* COMPAT 



coefficients of the objective function; 
basic variables; 

matrix of dimension M by M corresponding to 

the inverse of the basic matrix; 

set of indices corresponding to the basic 

variables ; 

dual variables; 

duration of the schedule; 

sequence of cargoes to pick up; 

index of the variable leaving the basis; 

index of the variable entering the basis; 

search direction; 

maximum feasible step size; 

ship number; 

objective function value; 
movement requirements ; 

travel distances from current ship ports to 
POEs (ports of embarkation) ; 
travel distances between POEs and PODs 
(ports of disembarkation) ; 
ship speed; 

number of days to complete the deployment; 
matrix of dimensions M by M that contains in- 
formation about the compatibility ship-cargo. 



* 

* 

* 

* 

* 

* 

* 

* 

* 

* 

* 

* 

* 

* 

* 

* 

* 

* 

* 

* 

* 

* 

* 

* 



* 



* 
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Subroutines : 



The subroutines and their objectives are: 

- SIMPLX - solves the revised simplex method; 

- RSTEP1 - step 1 of the revised simplex method; 

- RSTEP2 - step 2 of the revised simplex method; 

- RSTEP3 - step 3 of the revised simplex method; 

- PHIPRM - updates the "B" inverse matrix; 

- RINVRT - inverts the B matrix; 

- RDAYS - estimates an initial value for the number 

of days to complete the deployment; 

- SUBPR - generates feasible schedules; 

- RTIME - computes travel times. 

- RESULT - writes the output. 



Key parameters : 



NLOA 


- number 


of 


full 


shiploads of cargoes; 


NSH 


- number 


of 


ships 


t 




NPOE 


- number 


of 


ports 


of 


embarkation ; 


NPOD 


- number 


of 


ports 


of 


disembarkation ; 



Output : The output provides the following information 

- objective function value, 

- number of simplex iterations, 

- optimal (minimum) number of days to complete the 
deployment , 

- optimal primal solution, 



* - optimal dual solution, * 

* - ships' schedules, * 

* - sequence of cargos to pick up per ship, and * 

* - schedules' durations. * 

* * 

* Input / Output devices : * 

* Disk (MOVREQ DATA) input - device 07 * 

* Disk (TRAVEL DATA) input - device 08 * 

* Disk (FSTDST DATA) input - device 09 * 

* Disk (COMPAT DATA) input - device 11 * 

* Disk (SPD DATA) input - device 12 * 

* Disk (DEPLOUT DATA) output - device 10 * 

* * 

* * 

************************************************************ 

C MAIN PROGRAM 

* * 

* The master problem. * 

* * 

IMPLICIT DOUBLE PRECISION (A-H,0-Z), INTEGER (I-N) 
PARAMETER ( MM = 100, NN = 2000, KK = 2 ) 

DIMENSION B (MM) ,C(NN) ,SB(MM) ,U(MM) , BINV (MM , MM) ,IB(MM) , 
&WORK (MM) , XB (MM) ,XC0L(MM) ,MR(100,KK) ,SPD(MM) 

REAL A (MM , NN) ,TR(15,15) ,IT(30,15) 
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INTEGER MR , SHIP , XTIME (NN) , SEQ (NN , MM) , TAU , TAUL , TAUEST , 
&LSEQ (NN) , NLOA , NSH , NPOE , NPOD , CTSHIP 
LOGICAL COMP AT ( MM , MM ) 

CHARACTER* 13 MOVREQ , TRAVEL , FSTTIM 
COMMON /UNITS/ NIN, NOUT 



C 

C Initialize variables. 

C 

DATA NLOA, NSH, NPOE, NPOD / 3, 2, 5, 5 / 

DATA A , B , C , XB , BINV / 200000* 0 ., 100 *1 . DO , 2000 * 0 . DO , 
&100*1.0D0, 10000*0. 0D0 / 

DATA SB, WORK, U, IB /100 * 0 . DO , 100 * 0 . DO , 100 * 0 . DO , 100 * 0/ 
DATA SEQ, LSEQ ,SPD / 200000*0 , 2000*0 , 100*0 . 0D0 / 
DATA MR, TR, IT, XTIME /200*0 , 225*0 ., 450*0 ., 2000*0 / 

NIN = 2 
NOUT = 6 
JOUT = 0 
TAU = 1 
TAUL = 1 
TAUEST = 1 



C 

C Read the data from the data files. 
C 
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READ (07,*) ( (MR(I,J) , J=l,2) ,I=l,NLOA) 

READ (08,*) ( (TR(I, J) ,J=l,NPOD) ,I=l,NPOE) 
READ (09, *) ( (IT (I, J) , J=l,NPOE) , 1=1, NSH) 

READ (11, *) ( ( COMPAT ( I , J ) , J=1,NSH) r I=l,NLOA) 
READ ( 12 , * ) (SPD(I) ,1=1, NSH) 



C 

C Estimate an initial value for TAU. 
C 



CALL RDAYS (NLOA , NPOE , NPOD , NSH , TAU , TAUL , TAUEST , SPD , 
&TR, IT) 

MD = TAU 

WRITE(NOUT, 8000) NLOA , NSH , NPOE , NPOD , TAU 
WRITE (NOUT, 8030) TAUL, TAUEST 



C 

C Convert input to number of columns (M) and number of 
C rows (N) in the "A" matrix. 

C 

M = NLOA + NSH 
ITER = 0 
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5000 DO 10 I = 1 , MM 
B ( I ) = 1.0D0 

SB ( I ) = 0.0D0 
U(I) = O.ODO 
IB ( I ) = 0 
WORK (I) = O.ODO 
XB ( I ) = l.ODO 
XCOL(I) = O.ODO 
DO 20 J = 1 , NN 
C(J) = O.ODO 
XTIME(J) = 0 
A ( I , J ) = 0. 

SEQ ( J , I ) = 0 
LSEQ(J) = 0 
20 CONTINUE 

DO 30 K=1 , MM 

BINV ( I , K) = 0 . DO 
30 CONTINUE 

10 CONTINUE 

CTSHIP = 0 
SHIP = 0 
K = 0 

N = 2 *NLOA + NSH 
C 

C Generate input for the revised simplex method. 
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c 

DO 40 I = 1 , M 
IB ( I ) = I 
DO 50 J = 1 , M 

IF (I .EQ. J) BINV ( I , J ) 
50 CONTINUE 

40 CONTINUE 

DO 60 I = 1 , M-NSH 
C(I) = 1.0D0 
60 CONTINUE 

C 

C Generate artificial variables. 

C 

DO 70 J = 1 , M 
DO 80 K = 1,M 

IF (J .EQ. K) A ( J , K ) = 1 
80 CONTINUE 

70 CONTINUE 

C 

C Generate surplus variables. 

C 

DO 90 J = 1 , M 

DO 100 K = M+l , 2 *M-NSH 



= 1 . DO 
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IF ( ( J .EQ. (K-M) ) .AND. ( J . LE . (M-NSH) )) 

& A ( J , K ) = -1. 

100 CONTINUE 

90 CONTINUE 

SUM = 0 . DO 
DO 110 1=1, M 
U(I) = -C(I) 

SUM = SUM + C (IB (I) ) *XB (I) 

110 CONTINUE 
OBJ = SUM 
1000 CONTINUE 

C 

C Strategy to choose for which ship the next schedule 
C will be generated 
C 

SHIP = SHIP + 1 

IF ( SHIP .EQ. NSH + 1 ) SHIP = 1 



C 

C Generate columns as needed by the master problem. 

C 

CALL SUBPR (U , XCOL , TAU , M , N , NLOA , NPOE , NPOD , NSH , MR , TR , 
&IT , A , COMPAT , IB , XB , SHIP , XTIME , K , SPD , SEQ , LSEQ , CTSHIP ) 
DO 120 I = 1, M 
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SUM 



0. ODO 



DO 130 J = 1, M 

SUM = SUM + BINV (I, J) * A ( J , K) 

130 CONTINUE 

SB (I) = SUM 
120 CONTINUE 

C 

C Perforin the revised simplex method. 

C 

CALL SIMPLX ( A , B , C , XB , BINV , SB , U , WORK , IB , OBJ , N , M , JOUT , 
&K, ITER) 

IF (OBJ .LT. 10.0D-4) THEN 
IF ( TAU .E Q. MD) THEN 
MD = MD + 10 
TAU = TAUL 
GO TO 5000 
END IF 
NT = 1 

DO 140 I = 1 , M 

IF (XB ( I ) .GT. 1.0D-3 .AND. XB ( I ) .LT. .9 ) NT=0 
140 CONTINUE 

WRITE (NOUT, 8010) NT 
IF (NT .EQ. 1) THEN 
GO TO 1100 
ELSE 
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IF ( ITER .GT. 4000) STOP 



GO TO 5000 
END IF 
END IF 



GO TO 1000 



C 

C Write the results. 

C 

1100 WRITE (NOUT, 8020) TAU 

CALL RESULT ( JOUT , XB , U , C , A , IB , M , N , OBJ , ITER , SEQ , LSEQ , 
&XTIME , NSH , NLOA) 

8000 FORMAT (2 OX, ’PROGRAM OUTPUT ' , / , 20X , ' ============== ' , 

&// , 6X, 12, IX, 'SHIPLOADS' ,3X,I2,1X, 'SHIPS' ,3X,I2, 

&1X, ’ POES ’ , 3X, 12 , IX, ' PODS ' , / , 6X , ' INITIAL ESTIMATED’ , 

&' VALUE = TAU1 = ’,12,//) 

8010 FORMAT(6X, 'NT =’ , 12) 

8020 FORMAT (6X, ' *** FINAL (OPTIMAL) TAU = ’,12) 

8030 FORMAT ( 6X , ' * * TAUL = ',I2,2X,', TAUEST = ’,12) 

STOP 

END 

SUBROUTINE SIMPLX ( A , B , C , XB , BINV , SB , U , WORK , IB , OBJ , 
&N,M, JOUT, K, ITER) 

* * 
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This subroutine perforins the revised simplex method. 



* 



* 



IMPLICIT DOUBLE PRECISION (A-H f O-Z) f INTEGER ( I-N) 
PARAMETER ( MM = 100, NN = 2000 ) 

DIMENSION XB (MM) , B ( MM ) ,C(NN) , BINV (MM , MM) ,SB(MM) , 

&U (MM) , WORK (MM) ,IB(MM) 

REAL A ( MM , NN ) 

INTEGER ELL,XTIME (NN) 

COMMON /UNITS/ NIN, NOUT 

JOUT = 0 

200 CONTINUE 
C 

ITER = ITER + 1 

IF ( JOUT .EQ. 1 ) RETURN 

CALL RSTEP2 (XB , SB , SIGB , ELL , M , JOUT) 

IF ( JOUT .EQ. 2 ) RETURN 

CALL RSTEP3 ( XB , C , B , BINV , A , WORK , OBJ , IB , ELL , K,N,M, ITER ) 
IF (OBJ .LT. 10.0D-4) THEN 
NT = 1 

DO 10 I = 1 , M 

IF ( XB ( I ) .GT. 1.0D-3 .AND. XB(I) .LT. .90 ) NT=0 
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10 



CONTINUE 



IF (NT . EQ. 1) THEN 
ITER = ITER + 1 
RETURN 
END IF 
END IF 

IF ( MOD (ITER, 10) . EQ . 0 ) CALL RINVRT (BINV , A , IB , 

&WORK , M , N) 

CALL RSTEP1 ( A , C , SB , U , BINV , IB , N , M , K , JOUT) 

GO TO 200 
END 



SUBROUTINE RSTEP1 ( A , C , SB , U , BINV , IB , N , M , K , JOUT) 



* This subroutine performs the step one of the revised 

* simplex method 



IMPLICIT DOUBLE PRECISION (A-H,0-Z), INTEGER ( I-N) 
PARAMETER ( MM = 100, NN = 2000 ) 

DIMENSION C (NN) ,SB(MM) ,U(MM) , BINV (MM, MM) ,IB(MM) 
REAL A ( MM , NN ) 
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INTEGER XTIME(NN) 



COMMON /UNITS/ NIN, NOUT 
C 

TOLCON = l.D-6 
JOUT = 0 
C 

C Compute the duals . 

C 

DO 10 J=1 , M 
SUM = 0.D0 
DO 20 1=1, M 

SUM = SUM + BINV(I , J) *C (IB (I) ) 
20 CONTINUE 

10 U(J) = - SUM 
C 

K = 0 

VKMIN = 1.D30 
C 

DO 50 1=1, N 
C 

C Check if I is in IB. 

C 

DO 30 J=1 ,M 

IF ( I . EQ. IB(J) } GO TO 50 
30 CONTINUE 

SUM = C(I) 
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DO 40 J=1,M 



SUM = SUM + A(J, I) *U(J) 
40 CONTINUE 

IF (SUM .GE. VKMIN ) GO TO 50 
VKMIN = SUM 
K = I 
50 CONTINUE 

IF (VKMIN .LE. -TOLCON) GO TO 60 

JOUT = 1 

RETURN 

60 CONTINUE 

C 

C Form SB. 

C 

DO 80 1=1, M 
SUM = 0.D0 
DO 70 J=1 , M 

SUM = SUM + BINV (I,J) *A(J,K) 
70 CONTINUE 

80 SB (I) = SUM 

RETURN 
END 
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SUBROUTINE RSTEP2 (XB , SB , SIGB , ELL , M , JOUT) 



* * 

* This subroutine performs the step two of the revised * 

* simplex method * 



IMPLICIT DOUBLE PRECISION (A-H,0-Z), INTEGER ( I-N) 
PARAMETER ( MM = 100, NN = 2000 ) 

DIMENSION XB (MM) , SB (MM) 

INTEGER ELL 

COMMON /UNITS/ NIN, NOUT 



EPS = l.D-6 
ELL = 0 
SIGB = 1.D30 
DO 100 1=1, M 

IF ( SB ( I ) .LT. EPS) GO TO 100 
RATIO = XB (I) /SB (I) 

IF (RATIO .GE. SIGB) GO TO 100 
SIGB = RATIO 
ELL = I 

100 CONTINUE 

IF (ELL .EQ. 0) JOUT = 2 

RETURN 

END 
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SUBROUTINE RSTEP3 ( XB , C , B , BINV , A , WORK , OB J , IB , ELL , K , 
&N,M, ITER) 



* This subroutine perforins the step three of the revised 

* simplex method 



IMPLICIT DOUBLE PRECISION (A-H,0-Z), INTEGER ( I-N) 
PARAMETER ( MM = 100, NN = 2000) 

DIMENSION C (NN) ,XB(MM) , B ( MM ) , BINV (MM , MM) ,WORK(MM) 
REAL A (MM , NN) 

INTEGER ELL, IB (MM) 

COMMON /UNITS/ NIN, NOUT 

DO 10 1=1, M 
10 WORK ( I ) = A ( I , K) 

CALL PHIPRM(BINV, WORK, ELL, M) 

DO 30 1=1, M 

SUM = 0 . DO 
DO 20 J = 1 , M 

SUM = SUM + BINV (I , J) *B ( J) 

20 CONTINUE 

XB ( I ) = SUM 
30 CONTINUE 

IB (ELL) = K 
SUM = 0 .DO 
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DO 40 I = 1,M 

SUM = SUM + C (IB (I) ) *XB (I) 
40 CONTINUE 
OBJ = SUM 
RETURN 
END 



SUBROUTINE PHIPRM ( BINV , WORK , ELL , M) 



* This subroutine updates the BINV matrix. 



IMPLICIT DOUBLE PRECISION (A-H f O-Z), INTEGER (I-N) 
PARAMETER ( MM =100, NN = 2000 ) 

DIMENSION BINV(MM,MM) ,WORK(MM) 

INTEGER ELL 

COMMON /UNITS/ NIN, NOUT 

TOL = l.D-6 
SUM = 0.D0 
DO 10 I = 1 , M 

SUM = SUM + BINV (ELL, I) ‘WORK (I) 

10 CONTINUE 

YSUM = DABS (SUM) 

IF (YSUM .GE. TOL) GO TO 20 
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20 

30 

40 

50 

60 



WRITE (NOUT, 8000) SUM 



STOP 

CONTINUE 

SUM = 1. DO/SUM 

DO 30 I = 1 , M 

BINV ( ELL , I ) = SUM*BINV (ELL , I ) 

IF ( (BINV (ELL , I ) .LT. TOL) .AND. (BINV(ELL,I) 
& .GT. -TOL) ) BINV ( ELL , I ) = 0 . 0D0 
CONTINUE 
DO 60 J = 1 , M 

IF ( J .EQ. ELL) GO TO 60 
TEMP = 0 . DO 
DO 40 I = 1 , M 

TEMP = TEMP + BINV (J, I) *WORK(I) 

CONTINUE 

IF ( (TEMP .LT. TOL) .AND. (TEMP . GT . -TOL) ) 

& TEMP = 0.0D0 
DO 50 I = 1 , M 

BINV ( J , I ) = BINV ( J , I ) - TEMP* BINV (ELL , I ) 
IF ( (BINV ( J , I ) .LT. TOL) .AND. (BINV(J,I) 
& .GT. -TOL) ) BINV ( J , I ) = 0 . 0D0 

CONTINUE 
CONTINUE 
RETURN 
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8000 FORMAT (6X , ' **** ERROR **** NEW MATRIX WOULD BE', 
&* SINGULAR, INNER PRODUCT =’,G15.6) 

END 

SUBROUTINE RINVRT ( BINV , A , IB , WORK , M , N) 



* This subroutine reinverts the basis. 



IMPLICIT DOUBLE PRECISION (A-H,0-Z), INTEGER ( I-N) 
PARAMETER ( ZERO = 0.0D0, ONE = 1 . 0D0 ) 

PARAMETER ( MM =100, NN = 2000 ) 

DIMENSION BINV (MM , MM) , IB(MM), OMAT ( MM , MM ) , WORK ( MM) 
REAL A ( MM , NN ) 

COMMON /UNITS/ NIN, NOUT 

TOL = l.D-6 
DO 10 I = 1, M 

DO 20 J = 1, M 

BINV (I, J) = ZERO 
OMAT (I, J) = A ( I , IB ( J ) ) 

20 CONTINUE 

BINV (I, I) = ONE 
10 CONTINUE 

C 
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C Locate maximum magnitude element on or below the main 
C diagonal. 

C 

DO 30 K = 1, M 

IF ( K .LT. M) THEN 
IMAX = K 

AMAX = DABS (OMAT (K , K) ) 

KP1 = K+l 

DO 40 I = KP1 , M 

IF ( AMAX .LT. DABS (OMAT (I ,K) ) ) THEN 
IMAX = I 

AMAX = DABS (OMAT (I, K) ) 

ENDIF 

40 CONTINUE 

C 

C Interchange rows IMAX and K if IMAX is not equal to K 
C 

IF (IMAX .NE. K) THEN 
DO 50 J = 1, M 

ATMP = OMAT (IMAX, J) 

OMAT (IMAX, J) = OMAT ( K , J ) 

OMAT ( K , J ) = ATMP 
BTMP = BINV ( IMAX , J ) 

BINV ( IMAX , J ) = BINV ( K , J ) 

BINV ( K , J ) = BTMP 
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50 



CONTINUE 



ENDIF 

ENDIF 



C 

C Test for singular matrix. 

C 

IF (DABS (OMAT (K,K) ) .LT. 1.0D-6) THEN 
WRITE (NOUT, 8000) K, OMAT (K , K) 

ELSE 

DIV = OMAT (K , K) 

DO 60 J = 1, M 

OMAT ( K , J ) = OMAT(K, J) /DIV 

IF ( (OMAT (K , J ) .LT. TOL) .AND. (OMAT(K,J) 
& .GT. -TOL) ) OMAT ( K , J ) = 0.0D0 

BINV (K , J ) = BINV(K, J) /DIV 

IF ( (BINV (K , J ) .LT. TOL) .AND. (BINV(K.J) 
& .GT. -TOL) ) BINV (K , J ) = 0.0D0 

60 CONTINUE 

DO 70 I = 1, M 

AMULT = OMAT ( I , K) 

IF ( (AMULT .LT. TOL) .AND. (AMULT 
& .GT. -TOL) ) AMULT = 0 . 0D0 

IF ( I .NE. K) THEN 
DO 80 J = 1, M 

OMAT ( I , J ) = OMAT ( I , J ) - AMULT 
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& 



* OMAT (K , J) 



BINV ( I , J ) = BINV ( I , J ) - AMULT 

& * BINV (K, J) 

IF ( (BINV ( I , J ) .LT. TOL) .AND. 

& ( BINV ( I , J) .GT. -TOL ) ) 

& BINV ( I , J ) = 0.0D0 

80 CONTINUE 

ENDIF 

70 CONTINUE 

ENDIF 
30 CONTINUE 

8000 FORMAT ( ' * ERROR: BASIS IS S INGULAR ' , 14 , D1 5 . 6 ) 

RETURN 

END 

SUBROUTINE SUBPR (U , XCOL , TAU , M , N , NLOA , NPOE , NPOD , NSH , 
&MR , TR , IT , A , COMPAT , IB , XB , SHIP , XTIME , K , SPD , SEQ , LSEQ , 
&CTSHIP ) 



* This subroutine generates feasible (acceptable) columns * 



IMPLICIT DOUBLE PRECISION (A-H,0~Z), INTEGER (I-N) 
PARAMETER ( MM = 100, NN = 2000, KK = 2 , JJ = 2000 ) 
DIMENSION XCOL (MM) ,U(MM) ,UU(MM) ,V(MM) ,XB(MM) ,SA(MM) , 
&IB (MM) , SPD (MM) 
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REAL A ( MM , NN ) ,TR(15,15) ,IT(30,15) 

INTEGER VIND(MM) ,PRED(0:JJ) , LOAD(0:JJ) ,TIME(0:JJ) , 
& MR ( 100 , KK) , FROLD , TOLD , PATH ( 0 : MM ) , CURLD , COUNT , 

& LENGTH, LASTND, SHIP, TT , MLNGTH , XTIME ( NN ) ,CTSHIP, 

& CTSON ( 0 : MM) , LSEQ (NN) , TAU , STACK ( 0 : JJ ) ,TOP, 

& SEQ(NN,MM) , CTBACK ( 0 : MM) 

LOGICAL COMP AT ( MM , MM) 

DOUBLE PRECISION MIN , MINRC 
COMMON /UNITS/ NIN,NOUT 



C 

C Initialize the variables. 
C 

LIMIT = N 
1000 CTSHIP = 0 
2000 DO 10 I = 1 , MM 

SA ( I ) = 0.0D0 
UU(I) = 0 . DO 
V(I) = 0 . DO 
VIND(I) = 0 
10 CONTINUE 

NNEG = NLOA 
MINRC = 0.0D0 

C 

C Sort the duals. 
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c 

DO 20 I = l,NLOA 
UU(I) = U(I) 

20 CONTINUE 

IF ( N .GE. ( 2 * NLOA + NSH + 1) ) THEN 

DO 30 I = 1 , M 
DO 40 J = 1 , M 

IF ( (IB ( J ) .GT. ( 2 * NLOA + NSH)) .AND. 

& (A ( (NLOA + SHIP ) , IB ( J ) ) .NE. 1.D0) .AND. 

& ( XB ( J ) .GT. 0.5D0) ) THEN 

SA ( I ) = SA ( I ) + A ( I , IB ( J ) ) 

END IF 

40 CONTINUE 

IF ( S A ( I ) .LT. 1.0D0 ) UU(I) = -2.0D0 
30 CONTINUE 

END IF 
C 

C Check ship-cargo and ship-port compatibility. 

C 

DO 50 I = 1 , NLOA 

IF ( .NOT. COMP AT ( I , SHIP ) ) UU(I) = 99.0D0 

50 CONTINUE 

DO 60 I = 1 , NLOA 
MIN = 0.1D-6 
COUNT = 0 
IND = 0 
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DO 70 J 



1 , NLOA 



IF ( UU(J) .LT. MIN ) THEN 
MIN = UU(J) 

IND = J 
COUNT = 1 
END IF 

70 CONTINUE 

IF ( COUNT .EQ. 0 ) THEN 
NNEG =1-1 
GO TO 4000 
END IF 
V(I) = MIN 
VIND(I) = IND 
UU(IND) = 99.0D0 
60 CONTINUE 

4000 CONTINUE 

C 

C The Modified Depth First Search Algorithm. 

C 

C 

C Create all nodes out of the source and include them 
C in a stack. 
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c 



DO 80 I = 0,JJ - 1 
TIME (I) = 0 
STACK (I) = 0 
LOAD (I) = 0 
PRED(I) = 0 
80 CONTINUE 

DO 90 I = 1 , MM 

XCOL(I) = 0 . DO 
PATH ( 1-1 ) = 0 

CTBACK(I-l) = 0 
CTSON(I-l) = 0 
90 CONTINUE 

LENGTH = 0 

LASTND = 1 

TOP = 0 
CURLD = 0 
FROLD = 0 
TOLD = 0 

DO 100 I = NNEG ,1,-1 
LOAD (NNEG -1+2) 
PRED ( NNEG -1+2) 
STACK (NNEG - I + 1 
TOP = TOP + 1 
LASTND = LASTND + 



= VIND(I) 

= 1 

) = (NNEG -1+2) 
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100 



CONTINUE 



C 

C Ma 
c 

3000 



110 



120 



130 



140 



in loop to search for feasible schedules. 



CURLD = STACK (TOP) 

IF ( LENGTH . EQ . NNEG ) THEN 
RCOST = 0.0D0 
DO 110 I = 1, LENGTH 

XCOL (LOAD (PATH (I) ) ) = 1.0D0 
CONTINUE 

XCOL (NLOA + SHIP) = 1.0D0 
DO 120 I = 1 , M 

IF ( XCOL (I) .EQ. 1.0D0 ) RCOST = RCOST + U(I) 
CONTINUE 

IF ( RCOST .GT. -1.0D-4) RCOST = 0 . 0D0 
IF ( (RCOST .LT. 0.0D0) ) THEN 

N = N+l 

DO 130 I = 1 , M 

A ( I , N) = XCOL (I) 

CONTINUE 

XTIME(N) = TIME (PRED (CURLD) ) 

DO 140 J = 1, LENGTH 

SEQ ( N , J ) = LOAD ( PATH ( J ) ) 

CONTINUE 

LSEQ(N) = LENGTH 
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RETURN 



END IF 
END IF 

IF ( CURLD .EQ. 0 ) THEN 

IF ( MINRC .LT. O.ODO ) RETURN 
CTSHIP = CTSHIP + 1 
SHIP = SHIP + 1 
IF ( SHIP .GT. NSH) SHIP = 1 
IF ( CTSHIP .EQ. NSH ) THEN 
WRITE (NOUT, 8000) 

TAU = TAU + 1 
WRITE (NOUT, 8010) TAU 
GO TO 1000 
END IF 
GO TO 2000 
END IF 

TOP = TOP - 1 

IF ( N .GE. NN ) STOP ’RUN OUT OF SPACE’ 

IF ( PRED( CURLD) .EQ. 1) THEN 
LASTND = LASTND - LENGTH 
LENGTH = 0 

END IF 
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IF ( CTSON( LENGTH) . EQ . CTBACK ( LENGTH ) ) THEN 

DO 150 I = 1, LENGTH 

IF ( PRED(CURLD) . EQ . PATH (I) ) THEN 

LASTND = LASTND - LENGTH + I 
LENGTH = I 
GO TO 5000 
END IF 

150 CONTINUE 

END IF 



5000 PATH (LENGTH + 1) = CURLD 
C 

C Compute the travel time to pick up another cargo. 

C 

FROLD = LOAD (PRED (CURLD) ) 

TOLD = LOAD (CURLD) 

TT = 0 

CALL RTIME (NPOE , NPOD , MR , TR , IT , FROLD , TOLD , SHIP , NSH , 
StNLOA , TT , SPD) 

TIME (CURLD) = TIME ( PRED ( CURLD ) ) + TT 



C 

C Verify if it is feasible to pick up another cargo. 
C 
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IF ( TIME ( CURLD ) . LE . TAU ) THEN 



170 



160 



180 



190 



CTBACK( LENGTH) = CTBACK ( LENGTH ) + 1 
LENGTH = LENGTH + 1 
CTSON (LENGTH) = 0 

DO 160 I = NNEG ,1,-1 

DO 170 J = 1, LENGTH 

IF ( VIND(I) .EQ. LOAD (PATH (J) ) ) GO TO 160 

CONTINUE 

LASTND = LASTND + 1 
LOAD ( LASTND ) = VIND(I) 

PRED( LASTND) = CURLD 
TOP = TOP + 1 

CTSON (LENGTH) = CTSON ( LENGTH ) + 1 
CONTINUE 

DO 180 I = LASTND, (LASTND - CTSON ( LENGTH ) +1 ), -1 
STACK (TOP) = I 
TOP = TOP - 1 
CONTINUE 

TOP = TOP + CTSON (LENGTH) 

ELSE 

LASTND = LASTND - 1 
RCOST = 0.0D0 
DO 190 I = 1, LENGTH 

XCOL ( LOAD ( PATH ( I ) ) ) = 1 . 0D0 
CONTINUE 

XCOL (NLOA + SHIP) = 1 . 0D0 
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200 



& 



210 



220 



230 



DO 200 I = 1,M 

IF ( XCOL(I) .EQ. 1.0D0 ) RCOST = RCOST + U(I) 
CONTINUE 

IF ( RCOST .GT. -1.0D-4) RCOST = 0.0D0 
IF ( RCOST .LT. 0.0D0 .AND. LENGTH . GT . 

INT (NLOA/NSH ) -1 ) THEN 

IF ( CTBACK( LENGTH) .EQ. 0 ) THEN 
N = N + 1 
DO 210 I = 1 , M 

A ( I , N) = XCOL(I) 

CONTINUE 

XTIME(N) = TIME (PRED (CURLD) ) 

DO 220 J = 1, LENGTH 

SEQ (N , J ) = LOAD ( PATH ( J) ) 

CONTINUE 

LSEQ(N) = LENGTH 
IF ( RCOST .LT. MINRC ) THEN 
MINRC = RCOST 
K = N 
END IF 
END IF 

DO 230 I = 1, LENGTH 

XCOL ( LOAD ( PATH ( I ) ) ) = 0 . 0D0 
XCOL (NLOA + SHIP) = 0.0D0 
CONTINUE 

ELSE 
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DO 240 I 



1 , LENGTH 



XCOL ( LOAD ( PATH ( I ) ) ) = 0 . 0D0 
240 CONTINUE 

XCOL (NLOA + SHIP) = 0.0D0 
END IF 

CTBACK (LENGTH) = CTBACK ( LENGTH ) + 1 
GO TO 3000 
END IF 

IF ( N .GE. LIMIT + 1 ) RETURN 
CTBACK (LENGTH) = 0 
GO TO 3000 

8000 FORMAT (// ,6X, ’ TAU NOT FEASIBLE, INCREASE TAU ') 

8010 FORMAT (/,6X, 'NEW TAU = ',14) 

END 

SUBROUTINE RTIME (NPOE , NPOD , MR , TR , IT , FROLD , TOLD , SHIP , 
&NSH , NLOA , TT , SPD ) 



* This subroutine calculates travel times. 



IMPLICIT DOUBLE PRECISION A-H,0-Z ), INTEGER ( I-N ) 
PARAMETER ( KK = 2, MM = 100, NN = 2000 ) 

DIMENSION SPD (MM) 

REAL IT (30,15) , TR (15,15) 

INTEGER MR ( 100 , KK) , TT , TOLD , SHIP , FROLD 
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TT 



0 



C Calculating the travel time. 

IF (FROLD .EQ. 0) THEN 

TT = IDNINTt (IT (SHIP , MR ( TOLD, 1) ) + TR ( MR ( TOLD , 1 ) , 

& MR (TOLD , 2 ) ) ) / (24. * SPD(SHIP))) 

ELSE 

TT = IDNINT ( ( TR ( MR ( TOLD , 1 ) , MR ( FROLD , 2 ) ) + 

& TR ( MR ( TOLD , 1 ) , MR ( TOLD , 2 ) ) ) / (24. * SPD(SHIP))) 

END IF 
RETURN 
END 



SUBROUTINE RDAYS ( NLOA , NPOE , NPOD , NSH , TAU , TAUL , 
&TAUEST , SPD , TR , IT) 



* This subroutine calculates an initial estimate of the * 

* number of days to complete the deployment. * 

* * 

PARAMETER ( MM = 100 ) 

INTEGER NPOE , NPOD , TAU , TAUL , TAUEST , NLOA , NSH 
DIMENSION TR (15,15) , IT(30,15) 

REAL TRMIN , IT , TR , ITMIN , SPMAX 
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DOUBLE PRECISION SPD(MM) 



COMMON /UNITS/ NIN,NOUT 

TAU = 1 
TAUL = 1 
TAUEST = 1 
C 

C Calculate the minimum distance to travel. 

C 

TRMIN = 99999999. 

DO 100 1=1 , NPOE 

DO 200 J=1 , NPOD 

IF ( TR ( I , J ) .LT. TRMIN) TRMIN = TR ( I , J ) 

200 CONTINUE 

100 CONTINUE 

C 

C Calculate the minimum travel distance from the initial 
C ships' locations to the POEs . 

C 

ITMIN = 99999999. 

DO 300 I = 1 , NSH 

DO 400 J = 1 , NPOE 

IF ( IT ( I , J ) .LT. ITMIN) ITMIN = IT(I,J) 

400 CONTINUE 

300 CONTINUE 
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c 

C Calculate the maximum traveling ships' speeds. 
C 

SPMAX = -1. 

DO 600 J = 1 , NSH 

IF ( SPD ( J ) .GT. SPMAX) SPMAX = SPD(J) 

600 CONTINUE 

C 

C Compute the average number of trips per ship. 

C 

NTR = NLOA / REAL (NSH) 



C 

C Calculate an estimate for TAU. 

C 

TAUL = INT ( ( ( ( (NTR * 2) - 1) * TRMIN ) + ITMIN) / 

& ( SPMAX* 24 . ) ) 

TAUEST = INT (15. 5 + ( 0.8 * TAUL )) 

TAU = INT ( (0.7 * TAUL) + (0.3 * TAUEST) ) 

RETURN 

END 
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SUBROUTINE RESULT ( JOUT , XB , U , C , A , IB , M , N , OB J , ITER , 



&SEQ , LSEQ , XTIME , NSH , NLOA ) 



* This subroutine writes the solution to the output file * 

* — — * 

IMPLICIT DOUBLE PRECISION (A-H, O-Z), INTEGER ( I — N ) 
PARAMETER ( MM =100, NN = 2000, ZERO = 0.0D0 ) 

DIMENSION U (MM) , C (NN) ,XB (MM) , IB (MM) 

INTEGER SEQ (NN, MM) , LSEQ (NN) , XTIME (NN) 

REAL A (MM,NN) 

COMMON /UNITS/ NIN, NOUT 

IF (JOUT .GE. 2) GO TO 80 
WRITE (NOUT, 8000) OBJ 
WRITE (NOUT, 8005) ITER 
WRITE (NOUT, 8010) 

C 

C Is X(I) basic? 

C 

DO 30 1=1, N 

DO 10 J=1 , M 

INDEX = J 

IF ( IB ( J ) .EQ. I) GO TO 20 
10 CONTINUE 

GO TO 30 
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20 



CONTINUE 



WRITE (NOUT, 8020) I, XB (INDEX) 

30 CONTINUE 

WRITE (NOUT, 8030) 

WRITE (NOUT, 8040) (I,U(I) , I=l,NLOA) 

WRITE (NOUT, 8050) (I,U(I) ,I=(NLOA+l) ,M) 

DO 70 1=1, N 
DO 40 J=1 , M 

IF (IB ( J ) .EQ. I) THEN 

IF ( (XB ( J ) .GT. l.D-2) .AND. (I . GT . 

& 2*NLOA+NSH) ) THEN 

DO 90 L=NLOA+l , M 

IF ( A(L, IB ( J) ) .GT. .9) 

& WRITE (NOUT, 9010) (L-NLOA) 

90 CONTINUE 

WRITE(NOUT, 9020) (SEQ (IB ( J) , K) , 

6c K = 1 , LSEQ (IB ( J) ) ) 

WRITE (NOUT, 8090) IB ( J) ,XTIME (IB ( J) ) 
END IF 
GO TO 60 
END IF 

40 CONTINUE 

C 

C X(I) is non basic. 

C 
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SUM = C(I) 



DO 50 J=1 , M 

SUM = SUM + A ( J , I ) *U ( J ) 

50 CONTINUE 

GO TO 70 

C 

C X ( I ) is basic . 

C 

60 CONTINUE 

70 CONTINUE 

80 CONTINUE 

IF ( JOUT .EQ. 2) WRITE (NOUT, 8070) 

IF ( JOUT .EQ. 3) WRITE (NOUT, 8080) 

RETURN 

C 

8000 FORMAT (//, 6X, ’OPTIMAL OBJECTIVE FUNCTION VALUE IS’, 
&F12 . 5) 

8005 FORMAT (//, 6X, 'NUMBER OF ITERACTIONS = ’,15) 

8010 FORMAT (//,17X, 'OPTIMAL PRIMAL SOLUTION’,/) 

8020 FORMAT ( 18X , ' X ( ' , 13 , ' ) =’,F14.7) 

8030 FORMAT (//, 18X, 'OPTIMAL DUAL SOLUTION',/) 

8040 FORMAT (18X,'U(',I2,') =',F14.7) 

8050 FORMAT ( 18X , ’ V ( ' , 12 , ' ) =’,F14.7) 

8070 FORMAT (//,6X, 'PROBLEM IS UNBOUNDED FROM BELOW') 

8080 FORMAT (// ,6X, 'PROBLEM HAS NO FEASIBLE SOLUTION') 



74 



8090 FORMAT (6X, 'DURATION OF SCHEDULE ’,14,' IS : ’,13,/) 

9010 FORMAT (/, 6X, ' SCHEDULE FOR SHIP : ', 12 ,/, 6X ,’ CARGOES 

&TO CARRY: ’ ) 

9020 FORMAT (6X, 14) 

END 
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